# You're Making Us Look Bad: Can Partisan Embarrassment Dampen Partisanship and Polarization?
# Elizabeth C. Connors, Taylor N. Carlson, and Steven W. Webster
# Public Opinion Quarterly 
# Study 2 Replication Code 

# Install packages if needed: 
#install.packages("survey")
#install.packages("ltm")
#install.packages("MASS")
#install.packages("stargazer")

# Load Packages 
library(survey)
library(ltm)
library(MASS)
library(stargazer)

# Set working directory to where you have study2.csv stored 
setwd() 

# Read in data 
d <- read.csv("study2.csv")

# Set up survey design

design2 <- svydesign(ids = ~caseid, weights = ~weight, data = d) # Includes full sample
design <- svydesign(ids = ~caseid, weights = ~weight, data = d[d$pid7!=4 & d$pid7!=8,]) # Excludes independents 

# Sample size excluding pure independents
dim(d[d$pid7!=4 & d$pid7!=8,])

# Cronbach's Alpha for embarrassing scenarios
emb_frame <- data.frame(d[d$pid7!=4 & d$pid7!=8, c("scandal", "extreme", "moderate", "misinfo", "uninfo")])

cronbach.alpha(emb_frame)

# Descriptive statistics on embarrassment reported in text 
emb_scen_tab <- matrix(NA, nrow=5, ncol = 4)
emb_scen_tab[1,] <- prop.table(svytable(~scandal, design))
emb_scen_tab[2,] <- prop.table(svytable(~extreme, design))
emb_scen_tab[3,] <- prop.table(svytable(~moderate, design))
emb_scen_tab[4,] <- prop.table(svytable(~misinfo, design))
emb_scen_tab[5,] <- prop.table(svytable(~uninfo, design))

colnames(emb_scen_tab) <- c("Not at all embarrassing", "A little embarrassing", "Moderately embarrassing", "Very embarrassing")
rownames(emb_scen_tab) <- c("Sex Scandal", "Too Extreme", "Too Moderate", "Believe/Spread Misinformation", "Uninformed about Politics")

emb_scen_tab

## By Party 

## Republicans 
emb_scen_tab_rep <- matrix(NA, nrow=5, ncol = 4)
emb_scen_tab_rep[1,] <- prop.table(svytable(~scandal+rep, design), 2)[, 2]
emb_scen_tab_rep[2,] <- prop.table(svytable(~extreme + rep, design), 2)[,2]
emb_scen_tab_rep[3,] <- prop.table(svytable(~moderate + rep, design), 2)[,2]
emb_scen_tab_rep[4,] <- prop.table(svytable(~misinfo + rep, design), 2)[,2]
emb_scen_tab_rep[5,] <- prop.table(svytable(~uninfo + rep, design), 2)[,2]

colnames(emb_scen_tab_rep) <- c("Not at all embarrassing", "A little embarrassing", "Moderately embarrassing", "Very embarrassing")
rownames(emb_scen_tab_rep) <- c("Sex Scandal", "Too Extreme", "Too Moderate", "Believe/Spread Misinformation", "Uninformed about Politics")

emb_scen_tab_rep

## Democrats
emb_scen_tab_dem <- matrix(NA, nrow=5, ncol = 4)
emb_scen_tab_dem[1,] <- prop.table(svytable(~scandal+dem, design), 2)[, 2]
emb_scen_tab_dem[2,] <- prop.table(svytable(~extreme + dem, design), 2)[,2]
emb_scen_tab_dem[3,] <- prop.table(svytable(~moderate + dem, design), 2)[,2]
emb_scen_tab_dem[4,] <- prop.table(svytable(~misinfo + dem, design), 2)[,2]
emb_scen_tab_dem[5,] <- prop.table(svytable(~uninfo + dem, design), 2)[,2]

colnames(emb_scen_tab_dem) <- c("Not at all embarrassing", "A little embarrassing", "Moderately embarrassing", "Very embarrassing")
rownames(emb_scen_tab_dem) <- c("Sex Scandal", "Too Extreme", "Too Moderate", "Believe/Spread Misinformation", "Uninformed about Politics")

# --- Create Figure 2

# Add CIs to Democrat Plot
emb_scen_tab_dem_SE <- matrix(NA, nrow=5, ncol = 4)

## Scandal
scandal_n <- sum(svytable(~scandal+dem, design)[,2])

scandal_p1 <- prop.table(svytable(~scandal+dem, design), 2)[1, 2]
scandal_p2 <- prop.table(svytable(~scandal+dem, design), 2)[2, 2]
scandal_p3 <- prop.table(svytable(~scandal+dem, design), 2)[3, 2]
scandal_p4 <- prop.table(svytable(~scandal+dem, design), 2)[4, 2]

emb_scen_tab_dem_SE[1,] <- c(sqrt(scandal_p1*(1-scandal_p1)/scandal_n),
                             sqrt(scandal_p2*(1-scandal_p2)/scandal_n),
                             sqrt(scandal_p3*(1-scandal_p3)/scandal_n),
                             sqrt(scandal_p4*(1-scandal_p4)/scandal_n))

## Extreme 
extreme_n <- sum(svytable(~extreme+dem, design)[,2])

extreme_p1 <- prop.table(svytable(~extreme+dem, design), 2)[1, 2]
extreme_p2 <- prop.table(svytable(~extreme+dem, design), 2)[2, 2]
extreme_p3 <- prop.table(svytable(~extreme+dem, design), 2)[3, 2]
extreme_p4 <- prop.table(svytable(~extreme+dem, design), 2)[4, 2]

emb_scen_tab_dem_SE[2,] <- c(sqrt(extreme_p1*(1-extreme_p1)/extreme_n),
                             sqrt(extreme_p2*(1-extreme_p2)/extreme_n),
                             sqrt(extreme_p3*(1-extreme_p3)/extreme_n),
                             sqrt(extreme_p4*(1-extreme_p4)/extreme_n))

## Moderate
moderate_n <- sum(svytable(~moderate+dem, design)[,2])

moderate_p1 <- prop.table(svytable(~moderate+dem, design), 2)[1, 2]
moderate_p2 <- prop.table(svytable(~moderate+dem, design), 2)[2, 2]
moderate_p3 <- prop.table(svytable(~moderate+dem, design), 2)[3, 2]
moderate_p4 <- prop.table(svytable(~moderate+dem, design), 2)[4, 2]

emb_scen_tab_dem_SE[3,] <- c(sqrt(moderate_p1*(1-moderate_p1)/moderate_n),
                             sqrt(moderate_p2*(1-moderate_p2)/moderate_n),
                             sqrt(moderate_p3*(1-moderate_p3)/moderate_n),
                             sqrt(moderate_p4*(1-moderate_p4)/moderate_n))

## Misinfo 

misinfo_n <- sum(svytable(~misinfo+dem, design)[,2])

misinfo_p1 <- prop.table(svytable(~misinfo+dem, design), 2)[1, 2]
misinfo_p2 <- prop.table(svytable(~misinfo+dem, design), 2)[2, 2]
misinfo_p3 <- prop.table(svytable(~misinfo+dem, design), 2)[3, 2]
misinfo_p4 <- prop.table(svytable(~misinfo+dem, design), 2)[4, 2]

emb_scen_tab_dem_SE[4,] <- c(sqrt(misinfo_p1*(1-misinfo_p1)/misinfo_n),
                             sqrt(misinfo_p2*(1-misinfo_p2)/misinfo_n),
                             sqrt(misinfo_p3*(1-misinfo_p3)/misinfo_n),
                             sqrt(misinfo_p4*(1-misinfo_p4)/misinfo_n))

## Uninfo 

uninfo_n <- sum(svytable(~uninfo+dem, design)[,2])

uninfo_p1 <- prop.table(svytable(~uninfo+dem, design), 2)[1, 2]
uninfo_p2 <- prop.table(svytable(~uninfo+dem, design), 2)[2, 2]
uninfo_p3 <- prop.table(svytable(~uninfo+dem, design), 2)[3, 2]
uninfo_p4 <- prop.table(svytable(~uninfo+dem, design), 2)[4, 2]

emb_scen_tab_dem_SE[5,] <- c(sqrt(uninfo_p1*(1-uninfo_p1)/uninfo_n),
                             sqrt(uninfo_p2*(1-uninfo_p2)/uninfo_n),
                             sqrt(uninfo_p3*(1-uninfo_p3)/uninfo_n),
                             sqrt(uninfo_p4*(1-uninfo_p4)/uninfo_n))

emb_scen_tab_dem_SE

# Store Figure 2

png("figure2.png",width=8,height=10,units="in", res=600)
par(mfrow = c(2, 1))

barCenters <- barplot(t(emb_scen_tab_dem), beside = T,
                      main = "Democrat Embarrassment by Scenario",
                      ylim = c(0,1),
                      names.arg = c("Sex Scandal", "Too Extreme", "Too Moderate", "Believe/Spread\nMisinformation", "Uninformed\nAbout Politics"),
                      col = c("gray0", "gray25", "gray50", "gray75"),
                      ylab = "Proportion")
legend(0, 1, legend = c("Not at all embarrassing", "A little embarrassing", "Moderately embarrassing", "Very embarrassing"), fill = c("gray0", "gray25", "gray50", "gray75"), bty = "n")
segments(barCenters, t(emb_scen_tab_dem) + 1.96*t(emb_scen_tab_dem_SE), barCenters, t(emb_scen_tab_dem) - 1.96*t(emb_scen_tab_dem_SE))


# ADD CIs to Republican Plot 

emb_scen_tab_rep_SE <- matrix(NA, nrow=5, ncol = 4)

## Scandal
scandal_n <- sum(svytable(~scandal+rep, design)[,2])

scandal_p1 <- prop.table(svytable(~scandal+rep, design), 2)[1, 2]
scandal_p2 <- prop.table(svytable(~scandal+rep, design), 2)[2, 2]
scandal_p3 <- prop.table(svytable(~scandal+rep, design), 2)[3, 2]
scandal_p4 <- prop.table(svytable(~scandal+rep, design), 2)[4, 2]

emb_scen_tab_rep_SE[1,] <- c(sqrt(scandal_p1*(1-scandal_p1)/scandal_n),
                             sqrt(scandal_p2*(1-scandal_p2)/scandal_n),
                             sqrt(scandal_p3*(1-scandal_p3)/scandal_n),
                             sqrt(scandal_p4*(1-scandal_p4)/scandal_n))

## Extreme 
extreme_n <- sum(svytable(~extreme+rep, design)[,2])

extreme_p1 <- prop.table(svytable(~extreme+rep, design), 2)[1, 2]
extreme_p2 <- prop.table(svytable(~extreme+rep, design), 2)[2, 2]
extreme_p3 <- prop.table(svytable(~extreme+rep, design), 2)[3, 2]
extreme_p4 <- prop.table(svytable(~extreme+rep, design), 2)[4, 2]

emb_scen_tab_rep_SE[2,] <- c(sqrt(extreme_p1*(1-extreme_p1)/extreme_n),
                             sqrt(extreme_p2*(1-extreme_p2)/extreme_n),
                             sqrt(extreme_p3*(1-extreme_p3)/extreme_n),
                             sqrt(extreme_p4*(1-extreme_p4)/extreme_n))

## Moderate
moderate_n <- sum(svytable(~moderate+rep, design)[,2])

moderate_p1 <- prop.table(svytable(~moderate+rep, design), 2)[1, 2]
moderate_p2 <- prop.table(svytable(~moderate+rep, design), 2)[2, 2]
moderate_p3 <- prop.table(svytable(~moderate+rep, design), 2)[3, 2]
moderate_p4 <- prop.table(svytable(~moderate+rep, design), 2)[4, 2]

emb_scen_tab_rep_SE[3,] <- c(sqrt(moderate_p1*(1-moderate_p1)/moderate_n),
                             sqrt(moderate_p2*(1-moderate_p2)/moderate_n),
                             sqrt(moderate_p3*(1-moderate_p3)/moderate_n),
                             sqrt(moderate_p4*(1-moderate_p4)/moderate_n))

## Misinfo 

misinfo_n <- sum(svytable(~misinfo+rep, design)[,2])

misinfo_p1 <- prop.table(svytable(~misinfo+rep, design), 2)[1, 2]
misinfo_p2 <- prop.table(svytable(~misinfo+rep, design), 2)[2, 2]
misinfo_p3 <- prop.table(svytable(~misinfo+rep, design), 2)[3, 2]
misinfo_p4 <- prop.table(svytable(~misinfo+rep, design), 2)[4, 2]

emb_scen_tab_rep_SE[4,] <- c(sqrt(misinfo_p1*(1-misinfo_p1)/misinfo_n),
                             sqrt(misinfo_p2*(1-misinfo_p2)/misinfo_n),
                             sqrt(misinfo_p3*(1-misinfo_p3)/misinfo_n),
                             sqrt(misinfo_p4*(1-misinfo_p4)/misinfo_n))

## Uninfo 

uninfo_n <- sum(svytable(~uninfo+rep, design)[,2])

uninfo_p1 <- prop.table(svytable(~uninfo+rep, design), 2)[1, 2]
uninfo_p2 <- prop.table(svytable(~uninfo+rep, design), 2)[2, 2]
uninfo_p3 <- prop.table(svytable(~uninfo+rep, design), 2)[3, 2]
uninfo_p4 <- prop.table(svytable(~uninfo+rep, design), 2)[4, 2]

emb_scen_tab_rep_SE[5,] <- c(sqrt(uninfo_p1*(1-uninfo_p1)/uninfo_n),
                             sqrt(uninfo_p2*(1-uninfo_p2)/uninfo_n),
                             sqrt(uninfo_p3*(1-uninfo_p3)/uninfo_n),
                             sqrt(uninfo_p4*(1-uninfo_p4)/uninfo_n))

emb_scen_tab_rep_SE

barCenters <- barplot(t(emb_scen_tab_rep), beside = T,
                      main = "Republican Embarrassment by Scenario",
                      ylim = c(0,1),
                      names.arg = c("Sex Scandal", "Too Extreme", "Too Moderate", "Believe/Spread\nMisinformation", "Uninformed\nAbout Politics"),
                      col = c("gray0", "gray25", "gray50", "gray75"),
                      ylab = "Proportion")
legend(0, 1, legend = c("Not at all embarrassing", "A little embarrassing", "Moderately embarrassing", "Very embarrassing"), fill = c("gray0", "gray25", "gray50", "gray75"), bty = "n")
segments(barCenters, t(emb_scen_tab_rep) + 1.96*t(emb_scen_tab_rep_SE), barCenters, t(emb_scen_tab_rep) - 1.96*t(emb_scen_tab_rep_SE))

dev.off()

## Caption sample sizes 
table(d$rep)
table(d$dem)

# --- Table 2

## Scandal (Column 1)
scandal_age <- svyolr(as.factor(scandal) ~ age, design=design)
summary(scandal_age)
scandal_age_table <- coef(summary(scandal_age))
pval_scandal_age <- pnorm(abs(scandal_age_table[, "t value"]),lower.tail = FALSE)* 2
scandal_age_table_p <- cbind(scandal_age_table, "p.value" = round(pval_scandal_age,3))
scandal_age_table_p

scandal_woman <- svyolr(as.factor(scandal) ~ woman, design=design)
#summary(scandal_woman)
scandal_woman_table <- coef(summary(scandal_woman))
pval_scandal_woman <- pnorm(abs(scandal_woman_table[, "t value"]),lower.tail = FALSE)* 2
scandal_woman_table_p <- cbind(scandal_woman_table, "p.value" = round(pval_scandal_woman,3))
scandal_woman_table_p

scandal_white <- svyolr(as.factor(scandal) ~ white, design=design)
#summary(scandal_white)
scandal_white_table <- coef(summary(scandal_white))
pval_scandal_white <- pnorm(abs(scandal_white_table[, "t value"]),lower.tail = FALSE)* 2
scandal_white_table_p <- cbind(scandal_white_table, "p.value" = round(pval_scandal_white,3))
scandal_white_table_p

scandal_educ_r <- svyolr(as.factor(scandal) ~ educ_r, design=design)
#summary(scandal_educ_r)
scandal_educ_r_table <- coef(summary(scandal_educ_r))
pval_scandal_educ_r <- pnorm(abs(scandal_educ_r_table[, "t value"]),lower.tail = FALSE)* 2
scandal_educ_r_table_p <- cbind(scandal_educ_r_table, "p.value" = round(pval_scandal_educ_r,3))
scandal_educ_r_table_p

scandal_interest <- svyolr(as.factor(scandal) ~ interest, design=design)
#summary(scandal_interest)
scandal_interest_table <- coef(summary(scandal_interest))
pval_scandal_interest <- pnorm(abs(scandal_interest_table[, "t value"]),lower.tail = FALSE)* 2
scandal_interest_table_p <- cbind(scandal_interest_table, "p.value" = round(pval_scandal_interest,3))
scandal_interest_table_p

scandal_pid_strength <- svyolr(as.factor(scandal) ~ pid_strength, design=design)
#summary(scandal_pid_strength)
scandal_pid_strength_table <- coef(summary(scandal_pid_strength))
pval_scandal_pid_strength <- pnorm(abs(scandal_pid_strength_table[, "t value"]),lower.tail = FALSE)* 2
scandal_pid_strength_table_p <- cbind(scandal_pid_strength_table, "p.value" = round(pval_scandal_pid_strength,3))
scandal_pid_strength_table_p

scandal_dem <- svyolr(as.factor(scandal) ~ dem, design=design)
#summary(scandal_dem)
scandal_dem_table <- coef(summary(scandal_dem))
pval_scandal_dem <- pnorm(abs(scandal_dem_table[, "t value"]),lower.tail = FALSE)* 2
scandal_dem_table_p <- cbind(scandal_dem_table, "p.value" = round(pval_scandal_dem,3))
scandal_dem_table_p

scandal_ideology <- svyolr(as.factor(scandal) ~ ideology, design=design)
#summary(scandal_ideology)
scandal_ideology_table <- coef(summary(scandal_ideology))
pval_scandal_ideology <- pnorm(abs(scandal_ideology_table[, "t value"]),lower.tail = FALSE)* 2
scandal_ideology_table_p <- cbind(scandal_ideology_table, "p.value" = round(pval_scandal_ideology,3))
scandal_ideology_table_p

## Extreme (Column 2)

extreme_age <- svyolr(as.factor(extreme) ~ age, design=design)
#summary(extreme_age)
extreme_age_table <- coef(summary(extreme_age))
pval_extreme_age <- pnorm(abs(extreme_age_table[, "t value"]),lower.tail = FALSE)* 2
extreme_age_table_p <- cbind(extreme_age_table, "p.value" = round(pval_extreme_age,3))
extreme_age_table_p

extreme_woman <- svyolr(as.factor(extreme) ~ woman, design=design)
#summary(extreme_woman)
extreme_woman_table <- coef(summary(extreme_woman))
pval_extreme_woman <- pnorm(abs(extreme_woman_table[, "t value"]),lower.tail = FALSE)* 2
extreme_woman_table_p <- cbind(extreme_woman_table, "p.value" = round(pval_extreme_woman,3))
extreme_woman_table_p

extreme_white <- svyolr(as.factor(extreme) ~ white, design=design)
#summary(extreme_white)
extreme_white_table <- coef(summary(extreme_white))
pval_extreme_white <- pnorm(abs(extreme_white_table[, "t value"]),lower.tail = FALSE)* 2
extreme_white_table_p <- cbind(extreme_white_table, "p.value" = round(pval_extreme_white,3))
extreme_white_table_p

extreme_educ_r <- svyolr(as.factor(extreme) ~ educ_r, design=design)
#summary(extreme_educ_r)
extreme_educ_r_table <- coef(summary(extreme_educ_r))
pval_extreme_educ_r <- pnorm(abs(extreme_educ_r_table[, "t value"]),lower.tail = FALSE)* 2
extreme_educ_r_table_p <- cbind(extreme_educ_r_table, "p.value" = round(pval_extreme_educ_r,3))
extreme_educ_r_table_p

extreme_interest <- svyolr(as.factor(extreme) ~ interest, design=design)
#summary(extreme_interest)
extreme_interest_table <- coef(summary(extreme_interest))
pval_extreme_interest <- pnorm(abs(extreme_interest_table[, "t value"]),lower.tail = FALSE)* 2
extreme_interest_table_p <- cbind(extreme_interest_table, "p.value" = round(pval_extreme_interest,3))
extreme_interest_table_p

extreme_pid_strength <- svyolr(as.factor(extreme) ~ pid_strength, design=design)
#summary(extreme_pid_strength)
extreme_pid_strength_table <- coef(summary(extreme_pid_strength))
pval_extreme_pid_strength <- pnorm(abs(extreme_pid_strength_table[, "t value"]),lower.tail = FALSE)* 2
extreme_pid_strength_table_p <- cbind(extreme_pid_strength_table, "p.value" = round(pval_extreme_pid_strength,3))
extreme_pid_strength_table_p

extreme_dem <- svyolr(as.factor(extreme) ~ dem, design=design)
#summary(extreme_dem)
extreme_dem_table <- coef(summary(extreme_dem))
pval_extreme_dem <- pnorm(abs(extreme_dem_table[, "t value"]),lower.tail = FALSE)* 2
extreme_dem_table_p <- cbind(extreme_dem_table, "p.value" = round(pval_extreme_dem,3))
extreme_dem_table_p

extreme_ideology <- svyolr(as.factor(extreme) ~ ideology, design=design)
#summary(extreme_ideology)
extreme_ideology_table <- coef(summary(extreme_ideology))
pval_extreme_ideology <- pnorm(abs(extreme_ideology_table[, "t value"]),lower.tail = FALSE)* 2
extreme_ideology_table_p <- cbind(extreme_ideology_table, "p.value" = round(pval_extreme_ideology,3))
extreme_ideology_table_p

## Moderate (Column 3)

moderate_age <- svyolr(as.factor(moderate) ~ age, design=design)
#summary(moderate_age)
moderate_age_table <- coef(summary(moderate_age))
pval_moderate_age <- pnorm(abs(moderate_age_table[, "t value"]),lower.tail = FALSE)* 2
moderate_age_table_p <- cbind(moderate_age_table, "p.value" = round(pval_moderate_age,3))
moderate_age_table_p

moderate_woman <- svyolr(as.factor(moderate) ~ woman, design=design)
#summary(moderate_woman)
moderate_woman_table <- coef(summary(moderate_woman))
pval_moderate_woman <- pnorm(abs(moderate_woman_table[, "t value"]),lower.tail = FALSE)* 2
moderate_woman_table_p <- cbind(moderate_woman_table, "p.value" = round(pval_moderate_woman,3))
moderate_woman_table_p

moderate_white <- svyolr(as.factor(moderate) ~ white, design=design)
#summary(moderate_white)
moderate_white_table <- coef(summary(moderate_white))
pval_moderate_white <- pnorm(abs(moderate_white_table[, "t value"]),lower.tail = FALSE)* 2
moderate_white_table_p <- cbind(moderate_white_table, "p.value" = round(pval_moderate_white,3))
moderate_white_table_p

moderate_educ_r <- svyolr(as.factor(moderate) ~ educ_r, design=design)
#summary(moderate_educ_r)
moderate_educ_r_table <- coef(summary(moderate_educ_r))
pval_moderate_educ_r <- pnorm(abs(moderate_educ_r_table[, "t value"]),lower.tail = FALSE)* 2
moderate_educ_r_table_p <- cbind(moderate_educ_r_table, "p.value" = round(pval_moderate_educ_r,3))
moderate_educ_r_table_p

moderate_interest <- svyolr(as.factor(moderate) ~ interest, design=design)
#summary(moderate_interest)
moderate_interest_table <- coef(summary(moderate_interest))
pval_moderate_interest <- pnorm(abs(moderate_interest_table[, "t value"]),lower.tail = FALSE)* 2
moderate_interest_table_p <- cbind(moderate_interest_table, "p.value" = round(pval_moderate_interest,3))
moderate_interest_table_p

moderate_pid_strength <- svyolr(as.factor(moderate) ~ pid_strength, design=design)
#summary(moderate_pid_strength)
moderate_pid_strength_table <- coef(summary(moderate_pid_strength))
pval_moderate_pid_strength <- pnorm(abs(moderate_pid_strength_table[, "t value"]),lower.tail = FALSE)* 2
moderate_pid_strength_table_p <- cbind(moderate_pid_strength_table, "p.value" = round(pval_moderate_pid_strength,3))
moderate_pid_strength_table_p

moderate_dem <- svyolr(as.factor(moderate) ~ dem, design=design)
#summary(moderate_dem)
moderate_dem_table <- coef(summary(moderate_dem))
pval_moderate_dem <- pnorm(abs(moderate_dem_table[, "t value"]),lower.tail = FALSE)* 2
moderate_dem_table_p <- cbind(moderate_dem_table, "p.value" = round(pval_moderate_dem,3))
moderate_dem_table_p

moderate_ideology <- svyolr(as.factor(moderate) ~ ideology, design=design)
#summary(moderate_ideology)
moderate_ideology_table <- coef(summary(moderate_ideology))
pval_moderate_ideology <- pnorm(abs(moderate_ideology_table[, "t value"]),lower.tail = FALSE)* 2
moderate_ideology_table_p <- cbind(moderate_ideology_table, "p.value" = round(pval_moderate_ideology,3))
moderate_ideology_table_p

# Misinformation (Column 4)

misinfo_age <- svyolr(as.factor(misinfo) ~ age, design=design)
#summary(misinfo_age)
misinfo_age_table <- coef(summary(misinfo_age))
pval_misinfo_age <- pnorm(abs(misinfo_age_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_age_table_p <- cbind(misinfo_age_table, "p.value" = round(pval_misinfo_age,3))
misinfo_age_table_p

misinfo_woman <- svyolr(as.factor(misinfo) ~ woman, design=design)
#summary(misinfo_woman)
misinfo_woman_table <- coef(summary(misinfo_woman))
pval_misinfo_woman <- pnorm(abs(misinfo_woman_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_woman_table_p <- cbind(misinfo_woman_table, "p.value" = round(pval_misinfo_woman,3))
misinfo_woman_table_p

misinfo_white <- svyolr(as.factor(misinfo) ~ white, design=design)
#summary(misinfo_white)
misinfo_white_table <- coef(summary(misinfo_white))
pval_misinfo_white <- pnorm(abs(misinfo_white_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_white_table_p <- cbind(misinfo_white_table, "p.value" = round(pval_misinfo_white,3))
misinfo_white_table_p

misinfo_educ_r <- svyolr(as.factor(misinfo) ~ educ_r, design=design)
#summary(misinfo_educ_r)
misinfo_educ_r_table <- coef(summary(misinfo_educ_r))
pval_misinfo_educ_r <- pnorm(abs(misinfo_educ_r_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_educ_r_table_p <- cbind(misinfo_educ_r_table, "p.value" = round(pval_misinfo_educ_r,3))
misinfo_educ_r_table_p

misinfo_interest <- svyolr(as.factor(misinfo) ~ interest, design=design)
#summary(misinfo_interest)
misinfo_interest_table <- coef(summary(misinfo_interest))
pval_misinfo_interest <- pnorm(abs(misinfo_interest_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_interest_table_p <- cbind(misinfo_interest_table, "p.value" = round(pval_misinfo_interest,3))
misinfo_interest_table_p

misinfo_pid_strength <- svyolr(as.factor(misinfo) ~ pid_strength, design=design)
#summary(misinfo_pid_strength)
misinfo_pid_strength_table <- coef(summary(misinfo_pid_strength))
pval_misinfo_pid_strength <- pnorm(abs(misinfo_pid_strength_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_pid_strength_table_p <- cbind(misinfo_pid_strength_table, "p.value" = round(pval_misinfo_pid_strength,3))
misinfo_pid_strength_table_p

misinfo_dem <- svyolr(as.factor(misinfo) ~ dem, design=design)
#summary(misinfo_dem)
misinfo_dem_table <- coef(summary(misinfo_dem))
pval_misinfo_dem <- pnorm(abs(misinfo_dem_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_dem_table_p <- cbind(misinfo_dem_table, "p.value" = round(pval_misinfo_dem,3))
misinfo_dem_table_p

misinfo_ideology <- svyolr(as.factor(misinfo) ~ ideology, design=design)
#summary(misinfo_ideology)
misinfo_ideology_table <- coef(summary(misinfo_ideology))
pval_misinfo_ideology <- pnorm(abs(misinfo_ideology_table[, "t value"]),lower.tail = FALSE)* 2
misinfo_ideology_table_p <- cbind(misinfo_ideology_table, "p.value" = round(pval_misinfo_ideology,3))
misinfo_ideology_table_p

# Uninformed (Column 5)

uninfo_age <- svyolr(as.factor(uninfo) ~ age, design=design)
#summary(uninfo_age)
uninfo_age_table <- coef(summary(uninfo_age))
pval_uninfo_age <- pnorm(abs(uninfo_age_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_age_table_p <- cbind(uninfo_age_table, "p.value" = round(pval_uninfo_age,3))
uninfo_age_table_p

uninfo_woman <- svyolr(as.factor(uninfo) ~ woman, design=design)
#summary(uninfo_woman)
uninfo_woman_table <- coef(summary(uninfo_woman))
pval_uninfo_woman <- pnorm(abs(uninfo_woman_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_woman_table_p <- cbind(uninfo_woman_table, "p.value" = round(pval_uninfo_woman,3))
uninfo_woman_table_p

uninfo_white <- svyolr(as.factor(uninfo) ~ white, design=design)
#summary(uninfo_white)
uninfo_white_table <- coef(summary(uninfo_white))
pval_uninfo_white <- pnorm(abs(uninfo_white_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_white_table_p <- cbind(uninfo_white_table, "p.value" = round(pval_uninfo_white,3))
uninfo_white_table_p

uninfo_educ_r <- svyolr(as.factor(uninfo) ~ educ_r, design=design)
#summary(uninfo_educ_r)
uninfo_educ_r_table <- coef(summary(uninfo_educ_r))
pval_uninfo_educ_r <- pnorm(abs(uninfo_educ_r_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_educ_r_table_p <- cbind(uninfo_educ_r_table, "p.value" = round(pval_uninfo_educ_r,3))
uninfo_educ_r_table_p

uninfo_interest <- svyolr(as.factor(uninfo) ~ interest, design=design)
#summary(uninfo_interest)
uninfo_interest_table <- coef(summary(uninfo_interest))
pval_uninfo_interest <- pnorm(abs(uninfo_interest_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_interest_table_p <- cbind(uninfo_interest_table, "p.value" = round(pval_uninfo_interest,3))
uninfo_interest_table_p

uninfo_pid_strength <- svyolr(as.factor(uninfo) ~ pid_strength, design=design)
#summary(uninfo_pid_strength)
uninfo_pid_strength_table <- coef(summary(uninfo_pid_strength))
pval_uninfo_pid_strength <- pnorm(abs(uninfo_pid_strength_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_pid_strength_table_p <- cbind(uninfo_pid_strength_table, "p.value" = round(pval_uninfo_pid_strength,3))
uninfo_pid_strength_table_p

uninfo_dem <- svyolr(as.factor(uninfo) ~ dem, design=design)
#summary(uninfo_dem)
uninfo_dem_table <- coef(summary(uninfo_dem))
pval_uninfo_dem <- pnorm(abs(uninfo_dem_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_dem_table_p <- cbind(uninfo_dem_table, "p.value" = round(pval_uninfo_dem,3))
uninfo_dem_table_p

uninfo_ideology <- svyolr(as.factor(uninfo) ~ ideology, design=design)
#summary(uninfo_ideology)
uninfo_ideology_table <- coef(summary(uninfo_ideology))
pval_uninfo_ideology <- pnorm(abs(uninfo_ideology_table[, "t value"]),lower.tail = FALSE)* 2
uninfo_ideology_table_p <- cbind(uninfo_ideology_table, "p.value" = round(pval_uninfo_ideology,3))
uninfo_ideology_table_p

# --- Appendix Table B1
d$age4 <- NA
d$age4[d$age <= 29] <- "18-29"
d$age4[d$age > 29 & d$age <= 44] <- "30-44"
d$age4[d$age > 44 & d$age <= 64] <- "45-64"
d$age4[d$age > 64] <- "65 and Over"
table(d$age4[d$pid7 != 4 & d$pid7 != 8])

b1 <- matrix(NA, nrow = 5, ncol = 4)
b1[1,] <- prop.table(table(d$age4[d$pid7 != 4 & d$pid7 != 8], d$scandal[d$pid7 != 4 & d$pid7 != 8]), 1)[,4]
b1[2,] <- prop.table(table(d$age4[d$pid7 != 4 & d$pid7 != 8], d$extreme[d$pid7 != 4 & d$pid7 != 8]), 1)[,4] 
b1[3,] <- prop.table(table(d$age4[d$pid7 != 4 & d$pid7 != 8], d$moderate[d$pid7 != 4 & d$pid7 != 8]), 1)[,4]
b1[4,] <- prop.table(table(d$age4[d$pid7 != 4 & d$pid7 != 8], d$misinfo[d$pid7 != 4 & d$pid7 != 8]), 1)[,4] 
b1[5,] <- prop.table(table(d$age4[d$pid7 != 4 & d$pid7 != 8], d$uninfo[d$pid7 != 4 & d$pid7 != 8]), 1)[,4]

colnames(b1) <- c("18-29", "30-44", "45-64", "65+")
rownames(b1) <- c("Scandal", "Extreme", "Moderate", "Misinformed", "Uninformed")

b1 <- round(b1, digits = 2)
b1


# --- Appendix Table B2
scandal <- svyolr(as.factor(scandal) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
extreme <- svyolr(as.factor(extreme) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
moderate <- svyolr(as.factor(moderate) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
misinfo <- svyolr(as.factor(misinfo) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
uninfo <- svyolr(as.factor(uninfo) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)

summary(scandal)
summary(extreme)
summary(moderate)
summary(misinfo)
summary(uninfo)

# --- Appendix Table B3
scandal <- polr(as.factor(scandal) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)
extreme <- polr(as.factor(extreme) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)
moderate <- polr(as.factor(moderate) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)
misinfo <- polr(as.factor(misinfo) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)
uninfo <- polr(as.factor(uninfo) ~ age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)

summary(scandal)
summary(extreme)
summary(moderate)
summary(misinfo)
summary(uninfo)


# --- Appendix Table B4

# Weighted
  # See footnote 10 in MS text, coefficient on embarrassment 
persuade <- svyolr(as.factor(persuade) ~ emb + age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
#summary(persuade)
persuade_table <- coef(summary(persuade))
pval_persuade <- pnorm(abs(persuade_table[, "t value"]),lower.tail = FALSE)* 2
persuade_table_p <- cbind(persuade_table, "p.value" = round(pval_persuade,3))
persuade_table_p

public <- svyolr(as.factor(public) ~ emb + age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
#summary(public)
public_table <- coef(summary(public))
pval_public <- pnorm(abs(public_table[, "t value"]),lower.tail = FALSE)* 2
public_table_p <- cbind(public_table, "p.value" = round(pval_public,3))
public_table_p

private <- svyolr(as.factor(private) ~ emb + age + woman + white + educ_r + interest + pid_strength + dem  + ideology, design=design)
#summary(private)
private_table <- coef(summary(private))
pval_private <- pnorm(abs(private_table[, "t value"]),lower.tail = FALSE)* 2
private_table_p <- cbind(private_table, "p.value" = round(pval_private,3))
private_table_p

# With multiple testing correction 
persuade_pvals <- persuade_table_p[, "p.value"]
persuade_pvals_adj <- p.adjust(persuade_pvals, method = 'bonferroni')
persuade_pvals_adj

public_pvals <- public_table_p[, "p.value"]
public_pvals_adj <- p.adjust(public_pvals, method = 'bonferroni')
public_pvals_adj

private_pvals <- private_table_p[, "p.value"]
private_pvals_adj <- p.adjust(private_pvals, method = 'bonferroni')
private_pvals_adj

# --- Appendix Table B5 

# Unweighted 
# See footnote 11 in MS text, coefficient on embarrassment 

persuade <- polr(as.factor(persuade) ~ emb + age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)
public <- polr(as.factor(public) ~ emb + age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)
private <- polr(as.factor(private) ~ emb + age + woman + white + educ_r + interest + pid_strength + dem  + ideology, data = d[d$pid7!=4 & d$pid7!=8,], Hess = TRUE)

summary(persuade)
summary(public)
summary(private)
## With multiple testing correction
# See footnote 11 in MS text with coefficient on embarrassment

persuade_table <- coef(summary(persuade))
pval_persuade <- pnorm(abs(persuade_table[, "t value"]),lower.tail = FALSE)* 2
persuade_table_p <- cbind(persuade_table, "p.value" = round(pval_persuade,3))
persuade_table_p
persuade_pvals <- persuade_table_p[, "p.value"]
persuade_pvals_adj <- p.adjust(persuade_pvals, method = 'bonferroni')
persuade_pvals_adj

public_table <- coef(summary(public))
pval_public <- pnorm(abs(public_table[, "t value"]),lower.tail = FALSE)* 2
public_table_p <- cbind(public_table, "p.value" = round(pval_public,3))
public_table_p
public_pvals <- public_table_p[, "p.value"]
public_pvals_adj <- p.adjust(public_pvals, method = 'bonferroni')
public_pvals_adj

private_table <- coef(summary(private))
pval_private <- pnorm(abs(private_table[, "t value"]),lower.tail = FALSE)* 2
private_table_p <- cbind(private_table, "p.value" = round(pval_private,3))
private_table_p
private_pvals <- private_table_p[, "p.value"]
private_pvals_adj <- p.adjust(private_pvals, method = 'bonferroni')
private_pvals_adj